function [I,BW] = read_nd2(path)

%% Default Arguments
if nargin<1
    cd ~/Desktop/Lim/Dropbox/Yeast_Clustering/Code/
    path = '../../../All_Lim_Microscopy/Yeast/090312_tvs_yTS15-yTS16_w303.nd2';
end

%% Open file
[IN] = bfopen(path);
[num_channels,num_positions,num_time] = parse_meta(IN);
    
%% Load Images
I = cell(num_channels,num_positions);
h_bar = waitbar(0,'Loading Images...');

for c=1:num_channels
    for p=1:num_positions
        for t=1:num_time
            I{c,p}(:,:,t) = IN{p}{c+(t-1)*3,1};
            
            %Status
            waitbar(sum([c-1 p-1 t].*[num_time*[num_positions 1] 1])/prod([num_channels num_positions num_time]),h_bar)
        end
    end
end
close(h_bar)

BW = cell(size(I));
options = optimset(optimset('fminbnd'),'TolFun',1,'TolX',1);

%% Thresholding
h_bar = waitbar(0,'Identifying Masks...');
for c=1:num_channels
    for p=1:num_positions
        for t=1:num_time
            %Contrast
            It{c,p}(:,:,t) = uint16(round(mat2gray(I{c,p}(:,:,t))*(2^16-1)));
            %2D De-noising (Wiener Filter)
            It{c,p}(:,:,t) = wiener2(It{c,p}(:,:,t));
        end
            %Background Subtraction
            se = strel('disk',floor(sqrt(size(I{c,p},1))));
            It{c,p} = imtophat(It{c,p},se);
        for t = 1:num_time
            BWt = It{c,p}(:,:,t);
            %Otsu's method
            thresh = graythresh(BWt);
            BWt=mat2gray(BWt)>thresh;
            BWt=imfill(BWt,4,'holes');
            BWt=bwareaopen(BWt,10);
            %Optimization method
%             thresh = findThreshold(BWt,c,options);
%             BWt=BW{c,p}(:,:,t)>thresh;
%             BWt=bwareaopen(BWt,10);
%             if c>1
%                 c,p,t
%                 imshow(BWt)
%                 pause(.5)
%             end
%             BWt=bwmorph(BWt,'erode',1);
            
            BW{c,p}(:,:,t) = BWt;
            
            %Status
            waitbar(sum([c-1 p-1 t].*[num_time*[num_positions 1] 1])/prod([num_channels num_positions num_time]),h_bar)
        end
    end
end
close(h_bar)


%% Save image stacks
savedir = regexprep(path,'.nd2','/');
save([savedir 'I.mat'],'I');
save([savedir 'BW.mat'],'BW');

end